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The exact stationary state of an asymmetric exclusion process with fully parallel dynamics is 
obtained using the matrix product Ansatz. We give a simple derivation for the deterministic case by a 
physical interpretation of the dimension of the matrices. We prove the stationarity via a cancellation 
mechanism and by making use of an explicit representation of the matrix algebra we easily find closed 
expressions for the correlation functions in the general probabalistic case. Asymptotic expressions, 
obtained by making use of earlier results, allow us to derive the exact phase diagram. 

PACS number: 45.70.Vn. 



I. INTRODUCTION 

In this paper we describe the exact stationary state of a asymmetric exclusion process (ASEP) with fully parallel 
dynamics and open boundaries. This is a special case of the Nagel-Schreckenberg model for traffic flow Exact 
results have been known for some time for several update rules such as random sequential and sublattice parallel 
|^| H . For fully parallel dynamics and open boundary conditions mean field results have been obtained recently ||,[l0) ■ 
In a recent preprint, Evans et al. |p"i"|| presented an exact solution of this model using a site oriented matrix product 
Ansatz p2[ . Using an exlicit representation of the resulting algebra they calculated the current and density profile 
via generating function techniques. 

We will present a simple and physical derivation of the solution for deterministic bulk dynamics. This solution 
leads us to a bond oriented matrix product Ansatz resulting in a matrix algebra for stochastic bulk dynamics. Using 
an explicit representation of this algebra it is shown that difficult recursion relations can be circumvented and that 
integral expressions for the current and density profile can be given almost immediately. The resulting integrals can 
be calculated resulting in closed expressions similar to those of the random sequential case. 

The outline of the paper is as follows. In section || the model is defined and some notation is fixed. In section III 



we find and solve a simple recursion relation for the deterministic case. This solution can be naturally recast in the 



form of a matrix product through an interpretation of the dimension of the matrices, which is done in section IV. 
Section [v| is concerned with the formulation of a matrix product Ansatz for the general case of which solutions are 



presented in section VI . In the subsequent section it is shown that with a diagonalization pro cedur e closed expressions 



for the density profile and other correlation functions are obtained easily. Finally, in section VIII, the phase diagram 
is derived using asymptotic expression for the density profile and the current. 

II. THE MODEL 

The model is defined on a one dimensional lattice with L sites. Each site may be occupied with a particle or it may 
be empty. Configurations on the lattice are written as {r} = {n, . . . ,tl} numbering from left to right, where the 
presence or absence of a particle at site i is denoted by r, = 1 or Tj = 0. An asymmetric exclusion process (ASEP) 
is defined by imposing the following dynamics at each timestep t — ► t + 1 for the particles: If there is a particle at 
site i and if site i + 1 is empty, it hops to site i + 1 with probability p and remains at site i with probability 1 — p. If 
site i + 1 is occupied the particle at site i remains there with probability 1. This dynamics is applied to all particles 
at the same time, hence the name fully parallel. This dynamics is also known as the rule-184 cellular automaton 
which prescribes how the value of Tj at time t + 1 depends on the values of Tj_i, n and Tj+i at time t. Given the 
configuration {t} at time t, the configuration {t'} at time t + 1 is given by 

t[ = Pi-m-l(l - Tj) + (1 -pi)Ti(l - Ti+l) + TjTj+1, (i = 2, . . . ,L - 1), (1) 

where f>i are stochastic boolean variables with mean (pi) = p. The lattice is coupled to two reservoirs at the first and 
last sites. Particles may enter the system from the first reservoir if the first site is empty with rate a and they may 
leave the system with rate f3 into the second reservoir at the last site. Thus, 
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t[ =a(l-ri) + (l-pi)ri(l-T 2 )+r 1 r 2! (2) 
t' l =Pl-itl-i(1-t l ) + (1-$)t l , (3) 

where a and (3 are boolean variables such that (a) = a and = (3. Note that the dynamical rules have a manifest 
particle hole symmetry given by 

Ti -» l-Ti-i+l, 

Let us now associate to every configuration {r} a vector r defining an orthonormal basis of a Hilbert space. A 
state P(t) of the system will be any vector in this Hilbert space, i.e. 

P(i)=^P t (r 1 ,...,r i )r. (5) 
M 

The time evolution of such a state may be written as follows, 

P(f + l)=TP(t), (6) 

where T is called the transfer matrix. First of all one would like to know the stationary state of this model, i.e. the 
eigenstate of the transfer matrix with eigenvalue 1. As this state is time independent we will suppress from here on 
the temporal suffix of P. 

III. RECURRENCE RELATION 

For the deterministic bulk dynamics (p = 1) the bulk relations (|l|) simplify. Tilstra and Ernst pj have shown that 
in this case each configuration that can occur in the stationary state may be spatially divided into three parts, a free 
flow part, a jammed flow part and an interface of varying width. The free flow part is defined to be that part of the 
configuration up to the rightmost 00 pair and consists of isolated particles only. The jammed flow part starts with 
the leftmost 11 pair and contains isolated holes only. The jammed flow and the free flow can not overlap, but they 
may be separated by an interface consisting of a sequence of 10 pairs. Using this identification the dynamical rules 
become very simple. Denote the site of the last zero of the last 00 pair by / and the site of the first one of the first 
11 pair by j, then the rules are given by 



= n-i, (i = 2, ...,/), 

= n+u (i = j,---,L-l), (7) 

= n_i = 7- i+ i, (i = / + !,..., j-i). 



If there are no 00 pairs in a particular configuration we set / = 1 if n = 0, otherwise f — 0. Similarly, if there are no 
11 pairs, j = L if tl = L otherwise j = L + 1. In these cases the bulk relations (0) are not valid but one has to use 
the boundary relations (|^) and (|J). The equations of motion for deterministic bulk dynamics (p = 1) induces for 
the stationary state the equation 



P( n , . . . , Tf, (10)", T h ... t l ) = F( Ti \t 2 )J(t l \tl-i) x 

P(r 2 , ...,T f , (10)" +1 , Tj , . • • 7i_!) + J2 P (T2, ■ ■ ■ ,Tf, (10)«01(10r-«, T jt . . . T L ^) 

1 

Here, F and J are the transition rates for particles entering and leaving the system. They are given by 



(8) 



F(0|0) 


= 1 — a, 


J(l\l) 


= l-f3, 


F(1|0) 


= a, 


J(0|1) 


= P, 


F(0|1) 


= 1, 


J(1|Q) 


= 1, 


F(l|l) 


- o, 


J(0|0) 


= 0. 



(9) 



Iteration of equation <M) suggests the following Ansatz for the probabilities P, 
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P( n , ...,r f , (10)", T h ...,t l ) = i^(7i, . . . ,r / )Pi(n)P j (r J , . . .,r L ), (10) 

/-i 

P f (n,...,T/)=^ JJ F(T,|Ti+l), (11) 
i=l 

L-l 

i 5 j(Ti,..,r i )=y i - j+i n%+i|4 (12) 

where a;, y and P\{n) are to be determined and i?L is a normalization. Substituting this Ansatz into (j^) we find that 
P\{n) obeys the following recursion relation, 

n 

P x (n) - x-^Ptin + 1) + (1 - a)(l - /?) $>a: 2 ) I W) n - p , (13) 

p=0 

where we have set Pi(0) = 1. Explicit consideration of the equations of motion for P(0 ... 0) and P(0 . . . 01) determines 
the ratio between x and y. A convenient choice that fulfills this relation is x = 0, y = a. The recursion relation (pT 
can be solved easily and all the sums can be performed to give 

fl - a)8 n+1 - (1 - 0)a n+1 

P I (n) = (a/3)" U a>P 11 P,a . (14) 

p — a 

We thus have calculated the complete probability distribution function for the stationary state for deterministic bulk 
dynamics. The normalization in (|Io| ) is given by 

_ (l-a^-(l-^ 
p — a 

All expressions remain valid for a = (3 by taking the appropriate limits. 



IV. MATRIX PRODUCT 



We now will construct a matrix product representation of the stationary state. This will turn out to be useful for 
calculational reasons but it also helps us to find the solution for p < 1. First the vectors r are written as product 
states, 

L-l 

r = (g)v(r l ,r l+1 ). (16) 

i=l 

Next we will show how the form of the solution obtained in the previous section hints in the right direction. 

The exact form of Pi(n) as given by equation ( |T4"| ) suggests another view at the structure of the probability 
distribution function. This becomes clear by rewriting ([l4]) as 

n n 

P(n) = £ P (2k) - {aP) 1 ' 2 Po(2k - 1), (17) 

k=0 k=l 

where Pq is given by 

P (fc) = {xy)- l Pf{Tf, . . . 1 T f+k )P- 3 (T f+k+U . . .,Tj), 

= (a(3) n a n - k/2 (3 k/2 . (18) 

Pf and Pj are defined in equations ([ll]) and ( p^ ) . Thus equation ([Tt]) tells us that the weight of the interface of width 
In is equal to a sum over the positions of a separator. Configurations to the left of this separator are regarded as 
belonging to the free flow part and configurations to the right as belonging to the jammed part. There are different 
prefactors for even and odd positions of the separator. After equation ([17|) is substituted into equation ( |l0| ) one 
may even place the separator inside the free flow or the jammed flow, as the resulting additional terms vanish. This 
suggests that the stationary state may be written as a matrix product, 
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(W\ 



F S 
J 



L-l 



\v)/z L , 



(19) 



which is to be understood as a normal matrix multiplication but where the matrix elements are tensored. The matrices 
F and J govern the free and jammed flow respectively, while S is a matrix for the separator. Indeed, we find 



F = 
J = 
S = 



/ xF(0|0)v(00) xF(0|l)v(01) 
^F(l|0)v(10) 

yj(l|0)v(01)\ 

v yJ(0|l)v(10) yJ(l\l)v(ll)) 

a/3v(01) 
-(a/3) 2 v(10) 

(W\ = (l,l,0,a), (V\ = (/3, 0,1,1) 



/3(l-a)v(00) /3v(01)\ 
a^v(10) ) ' 

av(01) 
a/3v(10) a(l-/3)v(ll) 



(20) 

(21) 

(22) 
(23) 



In other words, the matrix operator in ( |l9| ) has two types of binary indices; the explicit ones referring to the type of 
flow, free or jammed, are contracted, and the internal ones, the occupation numbers Tj, are tensored. 



V. ALGEBRA 



Having found the stationary state for deterministic bulk dynamics (p = 1) and its representation as a matrix 
product, we now reverse the problem and show that the stationary state may be found by using the matrix product 
Ansatz (MPA) technique also for stochastic (p ^ 1) dynamics. In the following we will derive a matrix algebra from 
the MPA for the ASEP, i.e. for arbitrary p. We will show that a finite dimensional representation for this algebra can 
be found for p = 1 by using the solution found in section IV. 
As usual for the MPA, the matrix product 




will be written as follows, 



(24) 



where (A, B, C, D) is a vector on the basis v(00), v(01), v(10), v(ll) with matrix valued entries. It is to be understood 
that each entry of the tensor product is bra-ketted with ri {W| and I^Or^- In equation ( pi| ) the tensored indices are 
explicit, while the contracted indices are implicit, this in contrast to (|l9|). In the following we will show that the 
transfer matrix for the ASEP with fully parallel update rules may be written as T = TZT L ~ 1 C, for which the following 
mechanism ensures that (E4j) is a stationary state, 



Aw\c 




\v), 



(25) 



(26) 



In order to find C, 1Z and T we introduce the following probability distribution, 

P(ti, ■ ■ . ,Tj_i; ai,Ti+i . . . ,tl), 



(27) 



which corresponds to a partial updated sequence with ti, . . . , Tj_i at time t+1 and r^+i , . . . , tl at time t. The variable 
(Tj can attain three different values, say 0, 1, 2. The value corresponds to a hole both at time t and t + 1, the value 
1 to a particle at time t, and the value 2 to a particle moving into site i at time t+1. The introduction of such a 
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third state cr, = 2 was originally done by Rajewsky et al M and is necessary to incorporate correctly the fully parallel 
update rule. The probabilities ( ^7|) correspond with the matrix product Ansatz in the following way 



P(. . .,n-i;ai,T i+ i, ...) = n(W\ ■ ■■Y(T i -2,n-i)Y(T i _ 1 ,(T i )Y(Ti,T i+1 ) ■ ■ ■ \V) 1 



(28) 



where Tj = <Ji mod 2. This also shows why we need the fifth matrix X = Y(0, 2) for the intermediate states in ( |25| ) 
and (|26|). The equations, which define T, for the probabilities ( p7|) are explicitly given by 



(29) 



These equations immediately determine the matrix T. In a similar fashion the boundary operators C and 1Z can 
be calculated and they are given by 



P(.. 


.,0;0,.. 


) = P(. 


• ;0,o,...), 




P(.. 


.,0;1,.. 


) = P(- 


.;0,i,...), 




P(.. 


.,0;2,.. 


) = PP( 


..;1,0,...), 




P(.. 


.,1;0,.. 


) = (1- 


p)P(...;l,0,...) 


+ P(...;2,0, 


P(., 


.,1;1,.. 


) = P(- 


.;1,1,...) + J P(.. 


.;2,i,...), 


P(., 


.,1;2,.. 


) = o 







£ = 



/l - a 

' 1 - a 
a 
a 
V 





1-pO 
1 
P 0/ 



n = 



a 13 o o o> 

1-/3 1 

1/30 

v0 1 -/30/ 



(30) 



With these definitions, equations (|26|) and (|29j) imply the following algebra, 

AA = AA, AB = AB, AX=pBC, BC = (l-p)BC + XA, BD = BD + XB, 
CA = CA, CB = CB, DD = DD, CX = pDC, DC = (1 - p)DC. 



(31) 



There are some further relations related to the product form of t which forbids that products like v(ri_i , r^) (g) v(l — 
TjjTi+i) occur in any physical quantity, 



AD = AC* = BA = BB = BX = CC = CD = DA = DB = DX = 0. 
The boundary conditions from (|2^) with the explicit values of the matrices C and 1Z ( ^0| ) become 



(32) 



(W\A = {l-a) {W\A 

o(W\B = (l-a)o(W|S 

^W\C = a a (W\A+{\-p) x (W\C, 

i(W\D = a (W\B+ i{W\D 

o(W\X = Pl (W\C 



A\V)o = A|V) + /3B|V)i 

B\V)i = (l-P)B\V)i+X\V) 

C\V) = d\V) o +0D\V)i 

D\V)i = (1-/3)^1^)1 



(33) 



Any solution to (^l|), (|3^ ) and (|3^ ) thus automatically is a solution for the stationary state via (|25| ) and (|2^). It is 
however not obvious that the cancellation mechanism of (^5|) and ( p6| ) is ap prop riate for this problem. Indeed, we 
will see that for the case of p < 1 we will need a slightly weakened version of (|T]) and (^) . 

By observation of explicit solutions for small system sizes we also have inferred the following relations between the 
matrices, 



DCB 
BCB 
BCA 
DCA 

with the following boundary conditions, 



a(3 ((1 - p)CB + DD+ pa/3D) , 
a/3 (AB + BD +pa(3B) , 
a/3 ((1 - p)BC + AA+ pa(3A) , 
a(3(l - p) {DC + CA + paPC) , 



(34) 



o(W\A = p0(l-a) o (W\ B\V)i = pa\V) , 

o(W\B = p$i{W\ D\V)x - pa{l-0)\V) u 

x{W\CA = aPWW\A+{l-p)i{W\C) DC\V) = a/3(D|^)i + (1 - p)C\V) ), 

i{W\CB = a/3( (W\B + i(W\D) BC\V) = a(3(A\V) + B\V)i). 



(35) 
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We believe ( |34| ) and (|35|) to be true for arbitrary system sizes. It turns out that these relations are particularly 
convenient to obtain a representation. Using this representation in ( |3l| ) and ( |33|) to find the hatted matrices then 
gives an easy proof of the stationarity of the Ansatz (|24|) . 
The relations ( |32|) can be easily fulfilled by writing 



A = A® [ n n ) , B = B 



Q j, ^i^o 
c = c ® ' i o) ' = ® (o 1 



(36) 



and similarly for the hatted matrices 



A = A® n „ , B = S 



X = A". 

The boundary vectors are then written as 



1 0/ ' w \0 1 

1 




(38) 



|V>o = |Vo)®(j), 1^ = 1^)® (J 
o(W\ = (Wo | 0(1,0), i<W| = (Wi|®(0,l). 

VI. REPRESENTATIONS 

A. Representation for 1 — p = (1 — a)(l — /3) 

Rajewski ei aZ. || have already remarked that the product form (|24|) is exact for ordinary numbers instead of 
matrices on the line 1 — p = (1 — ct)(l — f3). Indeed, there exist a one-dimensional representation given by 

A=(3(l-a), B=l, C = af3, D = a(l - /3), 
A = (3{l-af, B=l-a, C = a/3(l - a), V = a, X = ap, (39) 

with 

(W |=/3, (Wi| = l, |V ) = 1, |Vi)=a. (40) 

B. Representation for p = 1 

In the case of deterministic dynamics in the bulk, p = 1, a two-dimensional representation of the subalgebra given 



by (31) for general values of a and (3 can be found. The matrices A, B,C and D can be easily read off from the 



solution in section [V and are given by 



(13(1 -a) 0\ n _ (fl 1 
~ 1 j ' ° ~ I a 



c , = (a/3 -a(3\ v _ (0 



. (41) 
a/3 J ' " ~ ^0 a(l - /3) J ' 

The representation for the hatted matrices can then be found using the algebra ( |3l| ) and with its boundary conditions 
(E33|) and is given by 



G 



A 
C 
X 



(3(1 - a 




) 2 




a(3(l-a) 











V 



((3(1 -a) I- a 
^ 

a(3 
a 



a(3 a(l - (3) 
a 



with 



(42) 



0\o| = (/.0h 0V,.| = (.;.li. : ( f 1 V |Vi) = r 



(43) 



C. Infinite dimensional representation 

In appendix ^ it is explained that there exists a basis {e„, /„} on which the matrices take the following form, 



V = B(l-p) = a(3q 



(q 1 ••■\ 

q 1 
00?1 

q 



V 



A = C = a(3q 



/q •••\ 
1?00 

1 q 

1? 
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V: 



(44) 



where q = y/l — p. This representation does not lead to divergent sums if 



a,/3 > 1- y/l-p. 



(45) 



A representation that is valid for all values of a and (3 can also be constructed, but this particular one will be useful 
for us in the sequel. First of all we would like to diagonalize E = A + B + C + D to facilitate further calculations. To 
find the eigenvalues and eigenvectors of E, we define the following vectors, 



|z))o = f>"e„, |z»i = 5> n /»- 



(46) 



It will also be convenient to define the following parameters, 

p-a p- (3 

a = , o=—— , 

aq pq 

so that the boundary vectors may be expressed as (see appendix |a|) , 

\V) = K \^\b)) , \V)! = K\b)) 
1-p 

where k is defined such that the normalization is given by 

(W\V)o=P, i(W\V) 1 =a. 



(47) 



(48) 



(49) 



In appendix ^| it is shown that there exist vectors \z;z 1 ))± that are linear combinations of the vectors defined in 
(Eft) which have the following properties, 



E\z-z 1 ))+ = k+(z)\z ] z- 1 )) + = a(3q[z+- z +q + q l ^\z;z + , 
E\z; z- 1 ))^ = A_|z; z- 1 ))^ = aPq (q - q- 1 ) \z; z' 1 ))-, 



(50) 
(51) 



By writing the boundary vectors as linear combinations of these eigenvectors, see (B5) and (B7), the normalization 
can be expressed as, 
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where 



Z L = (o(W\ + t (W\) E L - X (\V) + |VOi) 

-^-(A+(z) - A_)A + (z) L - 1 if (z, a)K(z, b), (52) 

s |=l 47T1Z 

g(*,c) = - z^A-r- v c = a,b, (53) 
(z c)(z c) 

and 

. l-p-(l-a)(l-/3) l-afe 

K = — = . 

ap(l — p) p 
Expression ( ]52] ) can be rewritten using the identities in appendix |c[ We then find 



(54) 



where 



fr(a)-y = S,( a )-fr(6) 
p(a — p) a — o 



S L (c) = -(^(c)-«»-i(c)(? 2 -l)) =-(i2i(c)+pa/3Jii_i(c)), (56) 
p P 



(c<i) jf dz A £ 



|*|=i 47riz 



n + 1 

^A+(z) L if(z,c)(z-z- 1 ). (58) 



Note that the integral representation ([58]) for Rl(c) is only valid for c < 1. Under this condition it can be calculated 
to give (|57|). Equation ( |57|) however is valid for all values of c, which may be checked explicitly for small system sizes 
or by using another infinite dimensional representation for c > 1. 

In deriving the phase diagram we will need the large L behaviour of Rl(c), Expression ( |57|) for Rl{c) is similar to 
that of the ASEP with random sequential update, its asymptotics can be calculated similarly Hh3]. By identifying 



the terms that have the largest contribution to the sums we find that n ~ aL and m ~ 2/(1 — c) for c < 1, m ~ y2aL 
for c = 1, m ~ (c— l)aL/c for c > 1, which gives, 

Jfc(c)«^(^)V(l)*^ forc<l, (59) 
^MD'^ fi»c=l, (60) 



(1 - c" 2 )A + (c) L forol, (61) 



where 



p 2 (l-a) 



— ' A+(l) = a/3(l + v ^^) 2 , A + (a)=a(3 ^j _ . (62) 



2 + 9 + 9 1 ac(p — a) 

A+(6) is obtained from A+(a) by interchanging a and /3. 

VII. EXPRESSIONS FOR THE CURRENT AND DENSITY 

Using the algebra it is easy to derive the following expression for the current J £ = p(ri(l — Tj+i))i, 

Jl = P^- UW\ + !(W\) E l CE L - 1 - 2 (|V) + |U)i) (63) 
= paj3(^(l-2J L - 1 )+paj3^(l-J L - 2 ) S j , (64) 
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from which we find by induction, 



J L =pa/3%^(l- J L - X ). (65) 



The density profile is much harder to find from the algebra and our strategy will be to express all correlation function 
in terms of the eigenvectors of E. In doing so, the correlation functions are easily expressed as integrals over the 
unit circle and can be calculated exactly by the residue theorem or asymptotically vi a th e saddle-point method. We 
first demonstrate this for the current. Calculating the action of C on \z; using (Bl) and re-expressing it in the 

eigenvectors \z;z))± we find, 

C\z; z- 1 )) + = a(3q [(1 + zq^z)), - (1 + z^q)^- 1 ))^ 

= ^ A+tz^A ^ 1>>+ " |Z; ■ (66) 



pa(3k I dz A 



Inserting this into (|63j) we find that 



Jl = f —A + (z)^- L K(z, a)K(z, b) 

J\z\=l 4?rlz 

{aR L - X {a)-bR L - X {b)), (67) 



(a - b)ZL 
which indeed fulfills (|65|). 

In order to find the density profile we now calculate the two point correlator (TiTi+i)^ which is given by 

(nn +1 ) L = -L ( (W\ + x{W\) E i ~ 1 DE L ~ i ~ 1 (|V0o + \V)i) . (68) 

This can be given in an integral representation using 

D\z- z- 1 ))+ = aP(l -p) [z(l + zq)\z))i - *T X (1 + ^ 1 «)l^ 1 »i] • (69) 



Re-expressing (B9f) as a linear combination of eigenvectors using (C4) it then follows that 



(nn + i)L = / ^-k + {wf-'K{ w ,a){i + w -\)w- n x 



' ]Z A + (z) L -' l - 1 K(z,b)z(l + zq)z r ' 



z\=i 2?riz 

= ^^- L y- R L -m-2{a)Rm{b) + b -*J L (70) 

a 2 f3 2 iiq 2 ^ aq + q 2 + 1 
= -z > Rm{a)RL- m -2{b) + 1 Jl- (71) 

Here we have made use of the fact that the product A + (u>) I ~~ 1 A + (;z) L ~ l ~ 1 can be written as a sum in two ways. 
Equation ( |70| ) is useful for studying the right boundary while ( f7l| ) is more suited for the left boundary. The density 
profile is given by 

(Ti)L = {TiTi +1 ) L + JL/p (72) 

The easiest way to analyse the density profile is by looking at its lattice derivative 

a 2 B 2 ka 2 

t L (i) = (n +1 ) L - (n) L = i L!^R l _ 1 (a)R L - t - 1 {b). (73) 

The value of the current and the asymptotic behaviour of iz,(i) will determine the phase diagram. 
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VIII. PHASE DIAGRAM 



A. The case p = 1 

The two dimensional representation (^7J) for this case allow a simple evaluation of the current and density. The 
current takes two different values corresponding to a low (LD) and a high density (HD) region. 

• Low density phase LD 

Here a < f3 and the current and density profile are given by 

1 + a 



{Ti)L = 



a / 1-/3 



1 + a V P 



1 + ~l^*- r,i , (75) 



where £ 1 = — log(a//3) and r = L — i. The density profile is flat except near the right boundary where it falls 
of exponentially from its maximum value (tl)l to the bulk value. 

• Transition line from LD to HD 

On this line a = f3. The current is still given by ( |74| ) but the density profile becomes linear, 

, > o. ( 1 — a i \ , 

(rdL=—- 1 + y ) • (76) 

1 + a \ a L J 

• High density phase HD This phase is characterized by a > (3 and the current and density can be obtained 
from those in the low density phase by the particle hole symmetry (^J). They are given by, 

^ TTTr (77) 

fa)i=l^(l-(l-<*)e- ,/£ ). (78) 

where = — log(/3/a). Thus the density profile is flat except near the left boundary where it increases 
exponentially from its minimum value {t\)l to its bulk value. 

B. General values of p 



The current (67) may take three different values depending on the parameters a and (3. These values correspond 
to a low density, a high density and a so called maximum current phase. The density profile in these phases will be 
calculated and will give rise to a further discrimination of phases within the low density and high density phase. 
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_ ^ LD„ 


MC 


LDi 

HDi 


\ 

\ 

\ 

\ HD„ 

\ 

\ 

\ 



a 



FIG. 1. Phase diagram in the a ~ (3 plane. Phase boundaries are at a,/3 = 1 — y/1 — p and at a = (3 where the transition is 
discontinuous in the density. On the dashed line, given by 1 — p = (1 — a)(l — /?), the mean field solution is exact. 

• Low density phase LDj 



This phase is chararcterized by the values o>6>loro;</3<l — ^1 — p. The current and bulk density p in 
this phase is given by 



J_ = 



pa(3 



a(p — a) 



p = 1 — J I a. 



A + (a)+pa(3 p—a 2 
We find an exponential decay of the density profile with a length scale 

£ = ia ~ ■< 

with 

e 1 — 

The slope of the density profile is given by (with r = L — i), 

t L (i) = {b- b- x )qj- (l - e" 1 ^) e-< r - 1 M 

_ (l-q)(p-2/3 + /?) / _ 17 A 

(p-a')(l-/?) V )° ■ 



(79) 



(80) 



(81) 



(82) 



Since b > 1 the slope of the density profile is positive and the density approaches its bulk value from above. 
Transition line from LDi to LDn 



On this line, where 6=lor/3 = l — \f\ — p, diverges while £ a remains finite. The bulk values are the same 
as in phase LDi, but the slope of the density profile now becomes 



t L U) = ^— ( 1 - c^A e-^-^r- 1 / 2 , 



2<?J_ 
P\fcrn 
a(p-a) (l-p) 1/4 

(p - o?) V5F(1 - VT^p) 



1 - e 



\ e -(^-l)/£a r -l/2_ 



(83) 
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Low density phase LDn 

This phase is determined by a > 1 < b or a < 1 — y/l — p < j3. Here the slope of the density profile still has a 
power law correction to the exponential but the power is different. 

tT (i) - qJ ~ A + (o)(A + (a)-A + (b)) / _ m \ /c 3/2 (M) 

The slope changes sign at the curve a = b^ 1 or 1 — p = (1 — a)(l — (3). This is the curve on which the mean 
field solution is exact and a one dimensional representation of our algebra exists. On this line the stationary 
state completely factorizes and the density profile is flat. 

Transition line from LDn to the maximal current phase MC 

On this transition line the current and bulk value of the density are as in the maximum current phase. The 
slope of the density profile on this transition line is given by 

t L (i) = Jl^l(i/L)-^r-^. (85) 

4V7T 

Near the right boundary the slope decays algebraically as r~ 3 / 2 . Near the left boundary the slope of the density 
profile decays algebraically with a power of 1/2 but the amplitude is of order l/L. Thus up to order 1/L 
corrections the density profile may be regarded as flat near the left boundary. 

The maximal current phase MC 

This part of the phase diagram is characterized bya,6<lora,/3>l — ^/l — p. In the maximal current phase 
the current attains it maximum value which is first reached on its phase boundaries. Its value and the bulk 
value of the density are given by 

*- - A + (iTt^ - 5" - '4 < 86 > 

The slope of the density profile in this phase is given by 

= -ii^^/VAr 372 - (87) 

Since its slope is negative the density approaches its bulk value of p = 1/2 from above as i -1 / 2 near the left 
boundary and from below as r -1 / 2 near the right boundary. 

The high density phases HDi and HDn 

The behaviour of the density profile in the high density phases and on their phase boundaries can be obtained 
from those of the low density phase by the particle hole symmetry, see equation (||), 

Tj-l — > 1 — T r , 

a «-» p. 



Coexistence line 



This line is characterized bya = 6>lora = /3<l — y/1 — p. The length £ a = remains finite but £ diverges. 
On this line one finds a linear profile with a positive slope, 

... p — 2a + a 2 

= V^r- (89) 
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J 



In the limit of small rates, i.e. a = pa and f3 = pf3 and p — ► 0, we recover the results fo r the A SEP with random 
sequential update jjjl^] . By taking p — > 1 the results reduce to those derived in subsection VIII A . Our results are in 
perfect agreement with those of 

In all phases and phase boundaries the current and bulk density p satisfy the following relation which defines the 
fundamental diagram 

v 7 ! - 4 W (1 - p)) . (90) 

Following Kolomeisky et al. |l6| , we may understand the phase diagram qualitatively by considering the domain wall 
dynamics. In this picture two characteristic velocities are important, the domain wall velocity and the collective 
velocity. The collective velocity is the drift of the center of mass of a momentary local perturbation of the stationary 
state and is related to the current by, 

n 

V coU = Q- p J{p)- (91) 

This velocity changes sign at p = 1/2 where the current takes its maximum value. For positive domain wall velocity 
((3 > a) and a < 1 — ^/l — p, an increase of the left boundary density leads to an increase of the bulk density since 
Vcoii > 0. This happens until the left boundary density equals 1/2, or a = 1 — ^/l — p. At this point V co ii changes sign 
and a perturbation will no longer spread into the bulk. The system is in the maximal current phase and a further 
increase of the left boundary density does not lead to an increase of the bulk density. For (3 < a the system does 
not enter the maximal current phase because of the negative domain wall velocity. The overfeeding however still 
occurs, which implies that further increase of the left boundary density beyond 1/2 does not lead to changes of the 
characteristic length scales in the high density phase. This is seen in the divergence of the length scale £ . 

The correlations for the ASEP with fully parallel dynamics are much stronger than for other dynamics. This 
becomes apparent when considering the relation between the length scales £ a ,b and the curents in the high and low 
density phase. The lengths £ a ,b can be written as, 

C 1 ^ 1 , (92) 



G X =-log h — T - / max ). (93) 



where 

^T 1 = -lnp-l 

1 - J J, 

This is in contrast to the random sequential and sublattice parallel dynamics |6[|7|,|16| where 

^ = -log(-^-). (94) 

\ "max / 

In the latter cases this relation could be obtained directly by considering the domain wall as a biased random walker. 
We have no simple argument for the fluctations in the domain wall position that leads to ( p3| ) in the case of fully 
parallel dynamics. 



IX. CONCLUSION 



We have presented a stationary state solution of an asymmetric simple exclusion process with fully parallel dynamics. 
In the case of deterministic bulk dynamics the solution, obtained directly from the master equations, has the form 
of a product over two-dimensional matrices. In contrast to the ASEP's with other dynamics, the matrices depend 
on two sites, instead of one. In the general case the stationary state can still be written as a product over matrices, 
but of infinite size. The stationarity of this product state can be proven by means of a cancellation mechanism which 
is a bit weaker than in other cases. We have calculated the exact phase diagram using an explicit representation of 
the matrix algebra. In this way we could via a diagonalization procedure derive expression for the current and the 
density profile with relative ease. 

The results are independent of, and agree with those of Evans et al. [ pd| which were obtained by means of a different 
Ansatz. They prove the strength and the flexibility of the matrix product Ansatz, though until recently the fully 
parallel dynamical models have resisted solution. Even when the resulting algebra is cubic (as in the present paper) or 
quartic (in |ll]]) a representation could be obtained. Of course, now that the formalism has been set up, many other 



13 



properties of the stationary state can be calculated. Instantaneous correlation functions are relatively straightforward. 
As our representation includes probability distributions involving consecutive time steps, it is to be expected that 
the present formalism is capable in principle of producing time dependent correlation functions. A more difficult test 
of the formalism is the calculation of the distribution of travelling times, for which it is necessary to follow a single 
particle as it flows through the system. 
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APPENDIX A: INFINITE DIMENSIONAL REPRESENTATION 



As an example for finding infinite dimensional matrices we explicitly construct the representation used in the main 
text. First we choose the following vectors as a basis, 

{e ,/o,ei,/i,e 2 ,/2,...}, (Al) 

where 



g n = (a/3q) n (A - a0(l - p)) n g , e„ = Ag n , /„ = Cg„ 



(A2) 



Here q = ^/l — p and we choose go such that 

Dfo = af3(l-p)f , Bf = af3e . 



(A3) 



We then find the action of the matrices A,B,C and D on these vectors from ([34]) and (|3J]). For example, besides 
(A3) we find for n > 1, 



Df n = a/3g(/„_i + qf n ), 
Bf n = af3q~ 1 (e n -i + qe n ), 



so that V and B are indeed given by (|f4|). On the basis (Al) the boundary vectors are given by 

±—^(Vo\ = (V 1 \= K {l,b,b 2 ,b 3 ,...), 

(Wo | =K(l,a,a 2 ,a 3 ,...) , 

(Wi| - -^(Wo\B = nq- 1 1, a, a 2 , a 3 

PP VP 



? • • ■ / i 



where 



and 



P-OL , P~ P 

a= , 5— > 

aq pq 



p(l-p-(l-q)(l-/3)) 



(A4) 
(A5) 



(A6) 
(A7) 
(A8) 

(A9) 

(A10) 



This representation does not lead to divergencies if a, b < 1 or a, [3 > 1 — y/1 — p. There are many possibilities in 
choosing the set of basis vectors. We could for example define go in a different way than was done in QA3|). The 
representation chosen here has some advantages which are exploited in the main text. It is however possible to choose 
a representation which is valid for all values of a and b (Evans et al. give an explicit example of such a representation 
p|). See Derrida et al. [Q for a similar discussion. 
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It turns out that for this representation we can find hatted matrices satisfying the relations on the first line of ( |3l| ) 
but not those on the second line. It is however possible to relax the conditions of (|3l|) a little in the following way. 
Every matrix will be pre-multiplied by another matrix. In particular C or D will be pre-multiplied by B or D (or 
\{W\ at the boundary). We therefore do not have to statisfy the relations on the second line of (|3l| ) identically but 
only up to a term that vanishes when acted on by B, D or i(W\. Since i(W\ — q))i =0 and B\ — q))i = D\ — q))i = 0, 
this is the case if this term is a matrix of which the columns are multiples of | — q))i. Thus, instead of the algebra 
obtained from pfl) , we find a solution of the algebra implied by 





B 
C 
D 

















• 











\V)r L = 



B 

C 

b 





( A 




B 




C 




\D 



\V), 



(All) 



(A12) 




( A \ 
B 

C 

b 



(A13) 



A solution to this algebra still automatically gives rise to a stationary state. We then find in addition to 

/ 1 (-q)-l (-q)-* (-q)-3 ■ ■ \ 
-q 1 (-q)- 1 (-q)- 2 
A = C + p*aP (-9? -9 1 HT 1 

(-qr (-9F -i 1 

■J 



C = af3 



V ! 

/l-2p -pH-q)- 1 

q i-2p 

q 

o o 



(A14) 



~p 2 (-q)- 2 
-pH-q)- 1 
1 - 2p 

q 



~p 2 (.-q)- 3 

~p 2 (-q) 2 

-p'i-q)- 1 
l-2p 



\ 



•7 



(A15) 



B = V - paP 



/ 1 

' -q 

{-qf 

\-qf 



V 



V = V + pa/3, 



(A16) 



X = pT> + pa(3e 



2A 



■J 

/ 1 
-q 
{-q) 2 
l-q) 3 

V ; 



(A17) 



•7 



APPENDIX B: EIGENVECTORS 



It follows from a direct calculation that the actions of the matrices on the vectors defined by (Bff) are given by 
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A\z)) = a/3q ({z' 1 + q)\z)) - z^eo) , 
B\z)) 1 = a(i(l + zq- x )\z))o, 
C\z)) = a0q((z- 1 +q)\z)) 1 -z- 1 f o ), 
D\z)) 1 = at0q(q + z)\z))i. 

Taking linear combinations of different vectors \z))i and \z))o such that terms with eo and fo drop out jH 
the following eigenvectors of E. 

\z; z- 1 ))* = z\z)) - z-^z- 1 )), + V± (z)\z)) 1 - ^(z-^z^h, 

r] + (z) = zq ■ , ri^(z) = -q. 

z + q 



The eigenvalues corresponding to these vectors follow easily from (Bl) and are, 

A+(z) = a/3q (z + z^ 1 + q + q^ 1 ) , A_ = a(3q (q - q^ 1 ) . 
From now on we take \z\ = 1. The following relations hold for a < 1 and b < 1 which are equivalent to a, f3 > 1- 

k (z- z- 1 )(A+(z) - A_) 



( {W\ + 1 {W\)\z;z-% 



p(3 (z — a)(z 1 — a) 



( (W\ + i(W\)\z;z- 1 ))-=0. 



The vectors |V)o and \V)\ can be expressed in the eigenvectors using (CI), from which we get, 

<1: {Z ~ Z ~ V) \z; 0>+ = K- 1 ^ (\V) + \V)i - /3 2 «| - q)U) 



i =1 Airiz (z - b){z~ l - b) 



1-/3 



The third term on the right hand side is a null vector, i.e. i(W| — q))\ = and E\ — q))i = 0, and does not 
calculations. 



(Bl) 

we find 

(B2) 
(B3) 

(B4) 

■>/r=p. 

(B5) 
(B6) 

(B7) 
enter the 



APPENDIX C: IDENTITIES 



The following identity which is frequently used throughout this paper can be conveniently calculated (or looked up 
in Jl3] ) by writing the denominator of the integrand as a sum of two geometric series and using the residue theorem, 



dz (z k — z k )(z — z 
=1 Aitiz [z — c)(z~ x — c) 

In a similar fashion the following integral can be calculated for c < 1, 

-l\L, 



c < 1. 



z|=i 27riz 



L—n L-\-n 

EE 



-(2 + z + z-') L z 



(z — c)(z 1 — c) 



2L 



2L 



L—n—k—1 



k=0 k=0 

for < n < L. Specialising to n = 1 and rewriting the terms in the sum we find that 

,-l\2fn I „ | „-l\jV ■ - N 



dz [z — z y 



\z\=l ^iz 



(2 + z + z~ 1 ) N ( C <i) ^ [2N-p\ p+1 , 
■)(z-i-c) ^ V N ) N+l [i + C) ■ 



Another important identity that we use to express vectors in terms of the eigenvectors of E is 

^ A + "f) g -A ( |z;z_1))+ \^ z ^-) + - (l + zqMz-'l) 

w=i ^ (^|*»x - ^lOh) (d + ^h«.| - (1 + ^x^-l) 
= / 1 -(l- g 2 )|-g)) 1 «- (Z - 1 |. 



(CI) 



(C2) 



(C3) 



(C4) 



1G 



This again can be simply evaluated using the residue theorem and the fact that q < 1. 
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